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We study the tunable quantum Hall effects in a non-Abelian honeycomb optical lattice which 
is a multi-Dirac-point system. We find that the quantum Hall effects present different features 
with the change in relative strengths of several perturbations. Namely, the quantum spin Hall 
effect can be induced by gauge-field-dressed next-nearest-neighbor hopping, which together with 
a Zeeman field can induce the quantum anomalous Hall effect characterized by different Chern 
numbers. Furthermore, we find that the edge states of the multi-Dirac-point system represent very 
different features for different boundary geometries, in contrast with the generic two-Dirac-point 
system. Our study extends the borders of the field of quantum Hall effects in a honeycomb optical 
lattice with multivalley degrees of freedom. 
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I. INTRODUCTION 

The honeycomb lattice, which is the brick that builds 
graphene, plays a significant role in promoting new 
physics. The most remarkable feature of the energy band 
of the honeycomb lattice system is that the low-energy 
excitations display a linear dispersion relation 1 and are, 
thus, described by massless Dirac fermions 2 . Further¬ 
more, massless Dirac fermions can be tuned into mas¬ 
sive Dirac fermions by tunable perturbations. More im¬ 
portantly, the transition from the gapless to the gapped 
phase can result in new physics, such as the quantum 
anomalous Hall (QAH) effect in the well-known Haldane 
graphene model 3 and quantum spin Hall (QSH) effect in 
the Kane-Mele graphene model 4 . The former model pio¬ 
neers a new route to realize non-zero integer Hall conduc¬ 
tance without Landau levels and the latter model estab¬ 
lishes the foundation to breed topological insulators 5,6 . 
Comparison between the two models indicates the power 
of the internal spin degrees of freedom to produce new 
physics in honeycomb lattice systems. 

On the other hand, the development of laser and ultra¬ 
cold atom techniques gives rise to various optical lattices 
in which ultracold Fermi atoms can be trapped to sim¬ 
ulate the phenomena in condensed-matter systems 7-13 , 
such as the topological Mott insulator 14 (MI) and su¬ 
perfluid (SF)-MI transition 15,16 . More significantly, in a 
laser field with a specific configuration, the trapped ul¬ 
tracold atoms can feel effective Abelian or non-Abelian 
gauge fields 17,18 . The couplings between the trapped ul¬ 
tracold atoms and the artificial gauge fields are equiva¬ 
lent to various interactions in cold-atom systems, such as 
spin-orbital couplings 19,20 . 

In this paper, we study a non-Abelian honeycomb opti¬ 
cal lattice which is a two-dimensional multi-Dirac-point 
system. Compared with the usual honeycomb lattice, 
in which only two independent Dirac points emerge, the 
non-Abelian honeycomb optical lattice here has eight 


Dirac points. Namely, the dimension of the internal val¬ 
ley degrees of freedom is extended from 2 to 8, and the ex¬ 
tended valley degrees of freedom are induced by the cou¬ 
plings between trapped atoms and the non-Abelian gauge 
fields. Likewise, we introduce three kinds of perturba¬ 
tions to open a gap of Dirac points independently. Ex¬ 
plicitly, two perturbations including staggered sublattice 
potentials and gauge-field-dressed next-nearest-neighbor 
(NNN) hopping are time-reversal invariant. The third 
is Zeeman splitting, which breaks the time-reversal sym¬ 
metry of the system. We find that staggered sublattice 
potentials give a trivial gapped phase with no gapless 
edge states in the ribbon geometry, and the gauge-field- 
dressed NNN hopping results in the QSH effect, with the 
spin direction lying in the ribbon plane. Interestingly, 
we find that the edge states connecting different internal 
valleys present different features even though the over¬ 
all Chern number is 0 in the presence of both staggered 
sublattice potentials and Zeeman splitting. Furthermore, 
it is shown that cooperations among the three kinds of 
perturbations can induce QAH effects characterized by 
different Chern numbers. Likewise, edge states show in¬ 
ternal structures for different types of edge boundaries. 
The results indicate that many-Dirac-point systems rep¬ 
resent more abundant physics in comparison with the 
simple two-Dirac-points systems, and the robustness of 
the bulk-boundary correspondence is global, not local. 

The paper is organized as follows. In Sec. II, we in¬ 
troduce the non-Abelian honeycomb lattice model with 
three kinds of perturbations. In Sec. Ill, we study 
the topological phase transition driven by three kinds 
of perturbations by investigating the change in the valley 
Chern numbers defined at Dirac points. Furthermore, we 
study the relations between the change in Chern numbers 
and the emergence of edge states under different geome¬ 
tries. In Sec. IV, we discuss some important relevant 
issues and give a brief summary. 


II. MODEL 


(b) 


We start with a non-Abelian honeycomb lattice model 
proposed in Ref. 21 and introduce three kinds of perturba¬ 
tions into this model to open a gap at each Dirac point, 
respectively. The tight-binding Hamiltonian takes the 
form, 


H = H 0 + H U ( 1 ) 

with 

H() = -t bjUijdi + H.C., ( 2 ) 

(*.i> 

Hi = -h J2 (apljOi + bpijbi) 

« i,3 » 

+ S E (aj a * - &]&») + (a}(T z ai + b\a z bij . 

i i 

( 3 ) 

Here t (G) is the hopping amplitude between the (next) 
nearest neighbor link (i,j) (((i,j))), a\=(a\^,a\^j with 

aj a (ai a ) denoting the creation (annihilation) operator of 
a fermionic atom with spin a (up or down) at A sublat¬ 
tice i (a similar definition applies for sublattice B). The 
unitary operator Uij (E/T) is associated with the link con¬ 
necting the (next) nearest-neighbor lattice points —»• 
17 21 . The unitary operators coupling fermionic atoms to 
non-Abelian fields along each hopping path from sublat¬ 
tice to sublattice B are taken as U\ = e iai(Tx , EG = 1, 
E /3 = e ia3(J y . The three nearest neighbor hopping pathes 

are di = a (^,-5), d 2 = a(0,1), d 3 = a (-^,-5). 
Similarly, the unitary operators along each hopping path 
between the NNN lattice points are chosen as U[ = U'J = 
e ijl(Tz , U\ = e il2(7z , = Ug = e^ 3<7 -, with NNN 

hopping paths bi = —b4 = a (-\/ 3 ,0), b 2 = —bs = 

“(■^T’f)’ b 3 = _b 6 = a ( _ ^F>|)- Here, «s=i ,3 or 
7 s=i, 2,3 is the gauge flux, and cr x , a y , and a z are the 
Pauli matrices for spin. S denotes an on-site energy for 
sublattice A and B, and h is the Zeeman splitting. For 
simplicity, we choose t = 1 as the energy unit and the 
distance a between the nearest sites as the length unit 
throughout this paper. 

In the momentum space, Hamiltonian (1) can be writ- 
ten in the basis ^ k = (a kt ,& kt , ^kj., &kj.) T : 

# = E^o(k)+iMk)]V> k . (4) 

k 


Ho (k) = 


0 

-tPo (k) 
0 

~tP 2 (k) 


-tPS (k) 0 

0 -tPi (k) 
-tP-L* (k) 0 

0 -tP Q (k) 


Hi (k) = - (tiP 3 (k) - h) (a z ®t 0 ) + 6 (uq 0 
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FIG. 1 : (Color online) (a) Band structures for the Hamil¬ 
tonian Ho with a s =i,3 — f • (b) Constant energy surfaces 

corresponding to the second band. The Dirac points for the 
middle two bands lie at the K s -i,...,6 points and G s =i,... ,6 
points in the first Brillouin zone (BZ). Here, we only mark 
Ki and K4 and choose 4 ^ as the unit of k x and k y . 


and 

P 0 (k) =e- ikd2 , 

Pi (k) = ie“* k ' dl + e _ik ' d3 , 

P 2 (k) =ie~ tk dl -e“ ik d3 , 
p?, (k) = 2 sin(k-bg). 

s=l,2,3 

Here, we set <a s =i,3 = 75=1,2,3 = f • r z is the Pauli matrix 
to span two sublattices. To and gq are 2x2 unit matri¬ 
ces. The energy spectrum for the Hamiltonian 4Eo(k) is 
shown in Fig. 1 (a), where one sees that the whole en¬ 
ergy spectrum displays a “particle-hole” symmetry with 
respect to the Fermi energy 6^=0 and any two adjacent 
touching bands exhibit Dirac-type energy spectra. The 
lower (upper) two bands touch at six M points in the 
first Brillouin zone. The middle two bands touch at six K 
points, only two of which are independent, with Kq 4 = 

^(±^, 0 ), K 2j3i 5 j6 = ^(±^,±l) and G-points 

with G 1)4 = (±f,o), G 2 ,3,5,6 = 7 (±ir>±|) in 

the first Brillouin zone [shown in Fig. 1(b)]. Here, we 
choose ^ as the unit of k x and k y in all of the figures in 
this paper. 


III. QUANTUM HALL EFFECTS INDUCED BY 
PERTURBATIONS 

To open a gap at each Dirac point, Hi (k) in Eq. (4) 
is turned on. The scheme for simulation of the gauge 
field with 7 s=i, 2,3 = \ is presented in the Appendix. It 
is straightforward to check that nonzero G, S, and h can 
independently open gaps at the eight inequivalent Dirac 
points. However, gaps from different perturbations may 
have different topological features, and the competitions 
or cooperations between these perturbations determine 
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TABLE I: Coefficients of the Dirac Hamiltonian at eight Dirac 
points. 
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the overall topological properties of the system. In order 
to explicitly show the connections between the topolog¬ 
ical feature and the amplitudes of perturbations, we ex¬ 
pand the Hamiltonian in Eqs. ( 5 ) and (6) around each 
Dirac point and obtain the relevant low-energy Hamilto¬ 
nian as follows: 


(k) = (k)rfc + 4 m) (k)rfo + 4 m) (k )Vz ( 7 ) 

Here, 77^, rj y and r] z are the Pauli matrices in the effective 
band basis, m labels the Dirac point, and k is measured 
from the Dirac points. We summarize three components 
(daT^(k), d^ m ^(k), of d^ m ^(k) for all the Dirac 

points in Table I. The valley Chern number of the lower 
band around Dirac point m can be calculated by 22 


C m = i S gn(4 x d( m )(k) x 3 * B d< TO >(k)) z sgn( 4 "*>), (8) 


from which, we obtain 


C Ki /g 4 = -bgn ^ - h - , (9) 

Cjq/Gx = T^sgn , (10) 

Cg 2 /g 5 =T^sgn ^(5-3ii , (11) 

Ccs/Ge = ±Isgn u + fi T j . (12) 


here, C m can take a series of values, and the bulk-edge 
correspondence can represent abundant features due to 
the extended internal valley degrees of freedom. Further¬ 
more, edge states can show different features for different 
boundaries. 

In order to present an explicit picture of the bulk¬ 
boundary correspondence in the multi-Dirac point sys¬ 
tem, we calculate the spectra of Hamiltonian (1) with 
zigzag (armchair) boundaries along the y(x) direction 
and periodic boundaries along the x(y) direction. The 
results are shown in Fig. 2 and 3 . For a ribbon with 
zigzag boundaries, by comparing the magnitudes of the 
gaps opened by £i, S, and h at each Dirac point, we can 
identify that the points (G3,Gs) in the bulk Brillouin 
zone are projected to the same momentum k x = — 
in Fig. 2. Similarly, the points (G2,Ge), (Gi,K 4), and 
(Ki, G4) are projected to the momenta k x = 2v ^ 7r , 

and 4x ^ 7r , respectively. The same conclusions can also be 
drawn from Fig. 1 (b). Thus, it is convenient to divide 
the eight Dirac points into two groups with each group 
involving four points. Explicitly, group I includes points 
{G3, G2, G5, Ge}, and group II includes points {Gi, 
K l5 K 4) G 4 }. See Fig. 2 (a) for details. Furthermore, 
each group can be divided into two subgroups and each 
subgroup includes two Dirac points. Namely, two sub¬ 
groups in group I are {G 3 ,G 2 } and{G 5 , G 6 }, while two 
subgroups in group II are {Gi,Ki} and{K4,G4}. The 
topological properties of Dirac points can be character¬ 
ized by the valley Chern numbers Cg { and Ckj • In order 
to guarantee the bulk-boundary correspondence, we de¬ 
fine the joint valley Chern numbers C32, C56, Cn, and 
C44, in which C32 = Cq 3 + Cg 2 - Likewise, C56, C11, and 
C44 have similar definitions. For a ribbon with armchair 
boundaries, with a similar strategy, we find that groups 
I and II are the same as in the case of a ribbon with 
zigzag boundaries. However, we find from Fig. 3 (a) that 
the four Dirac points in each group mix together. We 
show that double-pairing and quadruple-pairing of Dirac 
points under different boundary conditions strongly in¬ 
fluence the transport behaviors along different types of 
edges. The information on the edge states can be ex¬ 
tracted from the distribution of the number and spin 
density of the atoms along the edges. Define the number 
and spin density of the atoms as follows, 


Having obtained the expressions for the valley Chern 
numbers, we now discuss the influences of the perturba¬ 
tions on the topological properties of the system. Gen¬ 
erally speaking, the topological properties are the global 
features of the system, and some physical quantities have 
simple corresponding relations to the total Chern num¬ 
bers, such as the Hall conductance cr^y = ~^h^ 2 m Cm 
in the electron system, and the mass conductance oc 
Cm in the ultracold atom system, with SC m denot¬ 
ing the change in Chern number with the applied pertur¬ 
bations. For the two-Dirac-point system, C m can only 
take thevvalue 0 or d=l, and the bulk-boundary corre¬ 
spondence is simple. For the multi-Dirac-point system 


4 n) (t) = (. a\ai{i)cToa n (i)\g ), ( 13 ) 

4 ” 1 (*) = \(aW n (i)crTOin(i)\g) ■ ( 14 ) 

Here, n^\i) and denote the number and spin of 

atoms at site i in the nth edge state, a = a,b marks 
two sublattices and r = x,y,z labels three components 
of spin. a n {i) = [a n5t (i),a n4 (i)] T . 

Now, we discuss the evolution of the bulk-boundary 
correspondence by turning on different perturbations. 
For a ribbon with zigzag boundaries, starting with the 
nonzero S, the edge spectrum is shown in Fig. 2(a), from 
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FIG. 2: (Color online) Band structures of the Hamiltonian 
(1) in the ribbon geometry with zigzag edges and 40 sites in 
the y direction, corresponding to (a) (5, £i,/i) = (0.1, 0,0), 
(b) (0.1,0.05,0), (c) (0.1,0.2,0), (d) (0.1,0,0.3) Here, k x is 
measured in units of 2n/3. 



FIG. 3: (Color online) Band structures of Hamiltonian (1) 
in the ribbon geometry with armchair edges and 60 sites in 
the x direction, corresponding to (a) ((5, ti,h) — (0.1, 0, 0), (b) 
(0.1,0.05,0), (c) (0.1, 0.2,0), and (d) (0.1, 0,0.3). Here, k y is 
measured in units of 27r/3. 


which no gapless edge states emerge. Correspondingly, 
the joint valley Chern number C32 = C56 = C\\ — C44 = 
0 . The system lies in a trivial state with the total Chern 
number Ctot = C 32 + C56 + C\\ + C44 in the half-filling 
case. When t\ is turned on, one can observe from Eq. 
( 9 - 12 ) that only Cg 2 an d Cg 5 change signs when t\ G 
[5/3, (5]. Then, we have C32 = C' G 2 |t 1 ><5/3-C'G 2 |t 1 <<5/3 = 
+1, and C 56 = C GB \ tl>S /3 - C G5 \ tl <6/3 = - 1 - Mean- 
while, two pairs of gapless edge states emerge [Fig. 2 (b)], 
and each pair of edge states connects each pair of Dirac 


points {G 3 ,G 2 } and {G 5 ,G 6 }, respectively. To show 
the features of edge states, we plot the distribution of 
n£\i) and in Figs. 4(al)-4(a4), and we can de¬ 

termine the propagation direction of the nth edge state 
through the velocity v = dE n (k x ) /dk x . The opposite 
propagation directions correspond to the opposite signs 
of the valley Chern numbers C 32 and C$q. The consis¬ 
tency between the valley Chern number C 32 , C 56 and 
the gapless edge states satisfies the bulk-boundary cor¬ 
respondence. The transport pictures of the edge states 
are schematically illustrated in Fig. 4(c), where the edge 
states for C 32 = +1 move along the clockwise direction 
and the edge states for C 56 = — 1 move along the coun¬ 
terclockwise direction. Furthermore, from Figs. 4(a2)- 
4(a4), we find that s£} z {i) and exponentially de¬ 

cay as the distance to the boundary with a positive¬ 
negative oscillation, while s^ x ^ y (i) and s^ x / y {i) expo¬ 
nentially and monotonously decay as the distance to the 
boundary. The s^(i) have similar behaviors. Thus, 
we clarify that the spin polarizations of edge states al¬ 
most lie in the a-b plane, edge states with opposite spin 
polarizations propagate at opposite velocities at each 
boundary of ribbon (see Fig. 4(c)), and states with 
t\ G [5/3,5] are QSH states. Compared with the stan¬ 
dard QSH state proposed in the Kane-Mele model 4 which 
involves four bands, the number of effective bands in our 
model is two, but the increased internal valley degrees 
of freedom compensate the minimum four-band require¬ 
ment to induce QSH states. In other multi-Dirac-point 
models 23,24 , similar QSH states might exist if the gaps 
can been opened by some special perturbations. Of more 
practical significance, Ref . 25 has studied the transport 
properties of a silicene model with two kinds of Rashba 
spin-orbit couplings 26 , where the QSH phase has similar 
features to our model. When t\ is tuned to be larger 
than 5, we find that Cki , Cg 4 , Ck 2 , and Cg ± in group 
II change their signs, and the joint valley Chern number 
Cu = C 44 = 0. Even though two pairs of edge states 
emerge [see Fig. 2(c)], they can be adiabatically tuned 
into the bulk. The topological properties are similar to 
the case with t% G [5/3,5]. For a ribbon with armchair 
boundaries, the joint valley Chern numbers C 32 and C 56 
merge into a single one, C 3 + C 5 + C 2 + Cq, as shown 
in Fig. 3(a). They become indistinguishable, and the 
joint valley Chern numbers can not be well defined. Cor¬ 
respondingly, the edge states also mix together in Figs. 
4(bl)-4(b4). These can be called quasiedge states, be¬ 
cause they cannot induce the net spin current but only 
accumulate atoms at the armchair edges. We schemati¬ 
cally show the features of quasi-edge states in Fig. 4 (d). 
Actually, the different behaviors of the edge states be¬ 
tween zigzag and armchair boundaries can be understood 
as follows. Zigzag boundaries have an A-sublattice ter¬ 
minal at one boundary and a B-sublattice terminal at the 
other boundary, while armchair boundaries have both A- 
and B-sublattice terminals at both boundaries. Figures 
4(al)-4(a4) show that the chirality of the edge state can 


5 





(bl) 


V. 

0.0 


(d) 


N 

0.05 


(bJ) 


armchair 


(al) 


(b3) 

0.05 

t (b4) 

3 V. 

M 

-1 

4 

U3) 

0.1 

1 

4 

0.3 

(M) 1 

► . j~V* 

4 

(b3) 5° 3 (b4) 

N 

8 



1 - t 


1 

"" 1 

5 





2 

3 ^ 

i IrH 

% 


■ 4 

4 


4 

. 4 

8 


-0.05 

-0..1 

r 

03 

1 -0.1 

\ -0.1 


U2) 


o.s 




i 

N 

U 

4 0 .0 i 


(bl) 


(c) 

r 


ligzag 


-e- 


-0- 


r«‘ 

it 


6 

7 1 '* 


S 

N 

8 

S 

7 ** 

8 



FIG. 4: (Color online) (al)-(a4) [(bl)-(b4)] Distributions of 
the atom density n a and three spin components (s ax , s ay , s az ) 
along a zigzag [armchair] ribbon with 100 lattice sites for four 
edge states in Fig. 2(b) [Fig. 3(b)], labeled 1, 2, 3, and 4. 
Here, a = a,b denote two A and B sublattices, (c) Schematic 
of the movement of atoms in edge states shown in Fig. 2(b). 
The spin-s^ component is nearly suppressed, while spin-s x 
and -s y components have high weights, (d) Schematic of the 
accumulations of atoms in edge states corresponding to Fig. 
3(b). The spin polarization is strongly suppressed, and no 
net-current can be driven. 


only be well defined when the edge state is localized on 
one kind of sublattice. This is the physical reason for the 
difference between ribbons with zigzag and ribbons with 
armchair boundaries. These features indicate that the 
transport behaviors of edge states from different types of 
edges can be very different in multi-Dirac-point systems 
even though the globally structures are characterized by 
a single topological number. 

Now, we turn to the effect of Zeeman splitting. Sup¬ 
pose that the system lies in a trivial state initially with 
5 0 (Fig. 2 (a)). Turning to ft, when ft > y/ 3 S, 

Cki, Cg 3 , Cg 4 , and Cg 5 change their signs. For 
a ribbon with zigzag boundaries, straight forward cal¬ 
culations show that the joint valley Chern numbers 
C 3 2/56 = C Gs/!>\h>V 36 ~ C G 3 / !i \ h < 73 S = -1 and @ 44/11 = 
OWK i \h >\/35 -Cgu/k, \h < V3S = + 1 ? and there are four 
pairs of gapless edge states, shown in Fig. 2 (d). Further 
analysis identifies that the two pairs of edge states con¬ 
necting Dirac points in group I with C32 = C50 = — 1 
propagate along the counterclockwise direction and the 
edge states connecting Dirac points in group II with 
C44 = C\\ = +1 propagate along the clockwise direction, 
as shown in Figs. 5 (c)- 5 (d). In Figs. 5 (al)- 5 (b 4 ), we 
see that and si™y(i) are 0 while is nonzero 

for n = 1...8, and the relevant spin polarization direc¬ 
tions of the edge-state are also shown in Figs. 5 (c)- 5 (d). 
Since the total Chern number is equal to 0, one can not 
distinguish this topological phase by measuring the Hall 


FIG. 5: (Color online) (al)-(b4) Distributions of the atom 
density n a and three spin components (s a cc, s ay , s az ) along a 
zigzag ribbon with 100 lattice sites for eight edge states in 
Fig. 2(d), labeled 1 to 8, respectively. Here, a = a : b denote 
two A and B sublattices, (c) Schematic of the movement of 
atoms in edge states in Fig. 2(d) labeled 1 to 4. (d) Schematic 
of the movement of atoms in edge states in Fig. 2(d) labeled 
5 to 8. 


conductance. Figure 2 (d) shows that the edge states are 
separated into two groups according to valleys with re¬ 
spect to k x = Thus we call these quantum anoma¬ 
lous valley Hall states. If other perturbations or inter¬ 
valley scatterings are introduced to break the topological 
structures in group I or II, the quantum anomalous val¬ 
ley Hall states become measurable. For a ribbon with 
armchair boundaries, the joint valley Chern number de¬ 
fined in the zigzag case shows the same behaviors, i.e., 
C32 = C56 = — 1, and C44 = Cu = +1. Therefore, the 
quadruple-pairing of them does not induce destructive 
mixtures and the spectra of the edge states are shown 
in Fig. 3 (d); they show features similar to those in Fig. 
2(d), due to having the same joint valley Chern numbers 
C32 = C56 = — 1 , and C44 = Cn = +1. 

Having obtained the individual picture of nonzero G 
and ft, we turn to nonzero G and nonzero ft combined. 
Then the Hall conductance or mass conductance becomes 
parameter dependent and is determined by the total 
Chern number. We plot the phase diagram about t\jt 
and h/t with S = 0.1 in Fig. 6, from which shows that 
the Chern number can be tchanged from 0 to 2 by t\ 
and ft in this multi-Dirac-point system. To clearly see 
the correspondence between the Chern number and the 
edge states of zigzag and armchair boundaries, we start 
with the parameters (S, G, ft) = (0.1, 0.0, 0.0) initially, 
and we have C32 = C56 = C\\ = C44 = 0 . When ( 6 , 
G, ft) = ( 0 . 1 , 0 . 03 , 0 . 11 ), only Cg 5 changes its sign and 
we have C 5 6 = 0^1^+^ - C Gs \ 6>3ti+ ^ h = - 1 . 
Thus, the system is in the QAH state with C to t = C56 = 
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FIG. 6: (Color online) Phase diagram in the ti — h plane. 
Here, we set S = 0.1 1. 


— 1. A pair of edge states emerges [see Figs. 7 (al) and 
7 (bl) for the cases of zigzag boundaries and armchair 
boundaries]. When (S , G, h) = ( 0 . 1 , 0 . 03 , 0 . 18 ), then 
Cki and Gq 4 also change signs. We obtain QAH states 
with C44 = C11 = +1 and the Chern number of the 
whole system G = G 56 + C44 + C\\ = —1 + 1 + 1 = +1. 
The corresponding three edge states are shown in Figs. 
7 (a 2 ) and 7 (b 2 ); two pairs have positive chirality and 
one pair has negative chirality. Furthermore, when ( 6 , 
ti, h ) = ( 0 . 1 , 0 . 1 , 0 . 07 ), Cg 2 changes sign and we obtain 
QAH states with G = +2. In Fig. 7 (a 3 ), we see that the 
edge states corresponding to G56 = — 1 and G23 = 1 cou¬ 
ple with each other and open a gap to destroy nontrivial 
features, while the edge states corresponding to C\\ = 1 
and C44 = 1 show nontrivial features, and retain the over¬ 
all Chern number G = +2. Note that even the leftward 
two valleys in Fig. 3 (b) have Chern numbers G56 = — 1 
and C23 = 1, but no gap is opened between them in com¬ 
parison with the edge states in Fig. 7 (a 3 ). The reason 
is that the time-reversal symmetry is conserved in Fig. 
3 (b), and this symmetry protects the QSH state. On the 
other hand, the time-reversal symmetry is broken in Fig. 
7 (a 3 ), and the two edge states can couple with each other 
and a band gap is opened. In Fig. 7 (b 3 ), it is straight¬ 
forward to check that the two pairs of edge states at the 
quadruple-pairing Dirac points {G3, G2, G5, Ge} give 
the Chern number G = G3 + G5 + G2 + Cq = 0, and 
the overall Chern number G = +2 is protected by an¬ 
other quadruple-pairing of Dirac points, {G x , K 1? K 4 , 
G 4 }, with Chern number C = 1 + 1 = +2. Note that 
in Fig. 7 (b 3 ), the gapless edge states in the right val¬ 
ley are fake and not stable, and increasing Zeeman split¬ 
ting can break these fake gapless edge states. Thus, the 
bulk-boundary correspondence is globally robust, how- 




FIG. 7: (Color online) (al)-(a3) [(bl)-(b3)] Band struc¬ 

tures of Hamiltonian (1) for the zigzag [armchair] ribbon 
geometry with 40 [60] sites in the y [x] direction with pa¬ 
rameters (5,ti,h) = (0.1,0.03,0.11), (0.1,0.03,0.18), and 
(0.1,0.1,0.07), respectively. These three cases show the edge 
states of the QAH effect with the Chern number of the whole 
system C-— 1, +1, and +2. Here, k x , and k y are measured 
in units of 2tt/3. 


ever, the edge states represent different features for dif¬ 
ferent boundaries. 


IV. DISCUSSIONS AND SUMMARY 

In principle, the largest Chern number for this eight- 
Dirac-point system is 4 22 . However, the phase diagram 
in Fig. 6 shows that the largest Chern number is 2. For 
the present three kinds of perturbations, we cannot ob¬ 
tain a G = 4 phase. If we introduce more perturbations 
to increase the numbers of parameters in Eqs. (8)-(ll), 
we believe that the (7 = 4 phase can be achieved. For ex¬ 
ample, the Chern numbers of K\ % G4, G2, and Gq change 
signs, while the other four do not change signs. Thus, we 
can obtain the G = 4 phase. We focus on three kinds 
of perturbations in the present paper and leave the real¬ 
izations of the higher Chern number phases for further 
studies. 

In the electron system, the Hall conductance is well 
defined as the response coefficient of the electron current 
about the gradient of the electric potential, i.e., the elec¬ 
tric field intensity. Similarly, in the ultracold atom sys¬ 
tem, the Hall conductance corresponds to the response 
coefficient of the atom mass current about the gradient 
of the optical trapping potential. Thus, the non-zero 
Hall conductance gives the experimental signatures of the 
mass current and accumulations of atoms at the bound¬ 
ary of the optical lattice when the trapping potential is 
modulated. 
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In summary, we study quantum Hall effects in a non- 
Abelian honeycomb optical lattice which is a multi-Dirac- 
point system. The Hall conductance has different values 
with the change of relative strengths of several pertur¬ 
bations. We find that gauge-field-dressed NNN hopping 
can induce the QSH effect, and a Zeeman field can in¬ 
duce the so-called quantum anomalous valley Hall effect, 
which includes two copies of quantum Hall states with 
opposite Chern numbers and counter-propagating edge 
states. The coexistence of both perturbations gives the 
nonzero Hall conductance which is characterized by dif¬ 
ferent Chern numbers. Our study shows the power of 
extended internal valley degrees of freedom in driving 
abundant quantum Hall effects in a non-Abelian honey¬ 
comb optical lattice. 


V. ACKNOWLEDGMENTS 

G. Liu was supported by NSF of China (No. 11247011), 
This work was supported by the NSFC under grants Nos. 
11434015, 61227902, 61378017, 11374354, 11174360, 
11575051, 51201057, NKBRSFC under grants Nos. 
2011CB921502, 2012CB821305, SKLQOQOD under 
grants No. KF201403, SPRPCAS under grants No. 
XDB01020300 and Hebei NSF (Grant No. A 2012202022 , 
A2015208024). 


Appendix A: Simulating the gauge-field-dressed 
next-nearest-neighbor hopping 

In this section, we focus on how to simulate the gauge- 
field-dressed NNN hopping, i.e., the first term in Eq. (3), 
in our model in optical lattices. The first term in the Eq. 
(3) actually corresponds to a specific spin-orbit coupling 
if we take 75 = 1 , 2,3 = ±§. It is well known that in the 
Haldane model 3 , a periodic vector potential applied to 
a two-dimensional honeycomb lattice with zero net flux 
through each unit cell. The spacial magnetic flux den¬ 
sity does not change the nearest-neighbor hopping am¬ 
plitude, but it really causes the NNN hopping to have 
a chirality and breaks the time-reversal symmetry [see 
Fig. 8 (a)]. Subsequently, Kane and Mele extended this 
mechanism to the time-reversal-invariant system, with 
up-spin and down-spin electrons having opposite chiral¬ 
ity when they are hopping between NNN sites 4 . This 
kind of spin-dependent effective magnetic field is the in¬ 
trinsic spin-obit coupling in the Kane-Mele model. In 
Fig. 8 , arrows between NNN sites denote the directions 
of positive phase hopping for down-spin electrons. The 
parameters 75 = 1,3 = \ and 72 = —f correspond to the 
spin-dependent effective magnetic field in Fig. 8 (b) and 
7 s=i, 2,3 = § correspond to Fig. 8 (c). It is surprising 
that these two kinds of effective magnetic fields can give 
similar results for the system, and the latter corresponds 
to the results in the text. For comparison, we also show 
the Haldane model or Kane-Mele model in Fig. 8 (a). 



FIG. 8: (Color online) Illustration of the contours of the spin- 
dependent effective magnetic field. Arrows between the next- 
nearest-neighbor sites denote the directions of positive phase 
hopping for down-spin electrons, (a) The Haldane model; . 
(b) the parameters 75 = 1,3 = f, 72 = — f, and </> — 7 t/2; (c) 
the parameters 7 s=i, 2,3 — § and (f) — 2 tt/3. 

The spacial periodic magnetic field in the Haldane 
model can be simulated by use of the laser-induced- 
gauge-field method in an optical lattice 27 . With a similar 
method, we can design two other kinds of effective mag¬ 
netic field. To this end, we consider a cold-atomic system 
with each atom having an A-type level configuration 11 . 
The ground states |1) and |2) are coupled to the ex¬ 
cited state 1 3} through a spatially varying standing-wave 
laser field, with Rabi frequencies Q p =flsm0e~ lSl and 
H c =Heos 0e~' lS ‘ 2 , respectively. With the rotating-wave 
approximation, the laser-atom coupling Hamiltonian is 
given by 

n ( 0 0 n p \ 

Him = -- ( 0 0 (Al) 

2 \n; n* c -2A ) 

with eigenstates (dressing states) 

| Xl )= e - iSl cos 6 > |1) -e~ iS2 sin 9 12) (A 2 ) 

|X 2 )= cos^sinfle - *' 5 ' 1 |1) + cos ipcos6e~ lS2 |2) — sin ip |3) 

(A3) 

^^SHi^sm^e - * 51 |1) + sin^cos0e _ * S2 |2) + cos<^|3) 

(A4) 

and eigenvalues \ n =i,2,3=0,% (A^\/A 2 +0 2 ). Here, 
the single-photon detuning A = CJ 3 — uj\ — cj p , with 
CJ 3 , cji, and l d p the intrinsic frequency of atom states 
13), atom state 11), and laser f} p , respectively. In 
the new basis space |x)={|Xi) , |X 2 >, |X 3 >), the 

/V 2 /V /V 

primary atom Hamiltonian H=^ J rH int (r)+K (r) 
can be rewritten as (—iJW—A) 2 +V with 

M the atom mass and A and V the matrix with 

matrix elements (y n (r) | Vy m (r)) and 

Vn,m =An (r) ^ n ,m+(Xn (r)| V (r) \xm (r)), respectively. 
One can see that in the new basis the atom can 
be considered as moving in gauge potential A, 
which corresponds to an effective magnetic field 

B e //=(V x A) — | (A x A) 27,28 . 

We focus on the subspace spanned by the two 
lower eigenstates {|xi), IX 2 )}, which is redefined by 
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|Xt) = IXi) and |X^> = |X 2 > in the spin language. This 
gives an effective spin-1/2 system. For the large detun¬ 
ing (A^>H) case, both state |yq) and state \x\) are stable 
under atomic spontaneous emission from the initial excite 
state 13). Furthermore, we assume the adiabatic condi¬ 
tion, which requires that the off-diagonal elements of the 
matrices A are smaller than the eigenenergy differences 
|A m —A n | (m, n=l,2,3) of states \xm}- Under this adia¬ 
batic condition and in the basis space { |xt) ? |xi)}> the 
gauge potential A becomes diagonal and has the form 29 

A = ( A o’ A* ) < A5 > 


with 


A t = -A| = h (VSi cos 2 e+VS 2 sin 2 0) . (A 6 ) 

Here we neglect the correction to nearest-neighbor tun¬ 
neling induced by the change in potential V (r) thanks 
to the large detuning approximation. 

Consider the configuration of two opposite- 
traveling standing-wave laser beams 27,29 , which 
take Rabi frequencies Q P =Q, sin fay — 0) e~ lklX and 
Q C =Q cos fay — 0) e lklX . The effective gauge potential 
is generated as A^=— A±=hki cos (2 k^y — 20) e x . Here 
fci=fcsin 0 i, k 2 = k cos 0 \, with k the wave-vector number 
of the laser and 0\ the angle between the wave vector 
and the e y axis. In order to obtain the special spin-orbit 
coupling in our model, the value of /c 2 should be selected 
appropriately. When the wave vector of the laser 
beams satisfies = 2lir with l = 1,2,3, •••, the laser 
beams are commensurate with the optical lattice. Fur¬ 
thermore, we can obtain a series of Peierls phase factors 
which satisfy our model. For example, we take k^a = 2i r. 
The Peierls phase factors for the nearest neighbor 
hopping in Fig. 8 (b) are 0, and those for the NNN 
hopping in Fig. 8 (b) are (p^ 5 =—(pQ 2 =c^V^ki cos (20) and 
^i 3 =( ^ 5 i =( ^ 24 =( ^ 46 =d 5 with a=±l representing the up- 
and down-spin. Considering the C 3 rotational symmetry 
of the hopping phase factors in Fig. 8 (b), we can rotate 
the vector potential A by ±27r/3 to obtain the other two 
vector potentials. Therefore, the total effective vector 
potential and magnetic field can be written as 


— sin ( 47 n/ — 20 ) 
(A 8 ) 

The total accumulated phases for nearest-neighbor hop¬ 
ping are 0, and those for the NNN hopping along the 
arrowed directions in Fig. 8 (b) are 


Bgjj =—aAirhki 


67 r 

2 sin (27 xy + 20 ) cos ( — 


^13 — ^35 ~ ^51 — 7^24 ~ — ^46 ~ — ^62 ~ OiVSki COS (20) 

(A9) 

The contours of the magnetic field with 0 = 7 t /2 for 
down-spin are also plotted in Fig. 8 (b). 

Comparing Fig. 8 (b) with Fig. 8 (c), we find 
that the phase factors have different signs along 
the sites (1 —3) direction and the latter does 
not have the C 3 rotational symmetry. So, along 
this direction, we select the Rabi frequencies 

n p =nsin (~^k 2 x - \k 2 y - 0 e i (r iklX+ ^ klV ) 

and cos — \k 2 y — 0 e *( 2 k ^ x + 2 k ^v) 

and still set k^a = 27r. The laser field along the 
(1 —3) direction brings the non-zero phase factor 
(fis = — (^46 = — afak\a cos (20). Finally, the total 
effective vector potential and magnetic field for the Fig. 
8(c) can be written as 


A e ff = hki cos ( 2 k 2 y — 20 ) — sin fay + 20 ) sin ^ V3k2 %^ 
- %/?thk\ cos fay + 20 ) cos (^/Sk 2 X^ e y (A 10 ) 


and 


~B^ff=a4:hki7r 


2 tt A . ( 2^3, 


2 cos ( —y + 20 ) sin 


+2 sin ( — y — 20 


(All) 


The phase factors along the arrowed direction in Fig. 8 (c) 
are 


A_ e j j‘ — (y,kk\ 


^cos (2k 2 y — 20 ) — cos (. k 2 y + 20 ) ( 


-Vshh sin fay + 20 ) sin (^/3k2X^j < 


/ r \ \^35 — ^51 — ^13 — ^46 — ^24 ~ 2 — OlVSki COS (20) 

yV3k 2 xjje x (A 12 ) 

The contours of the magnetic field with 0 = 27 t/ 3 for 
(A?) down-spin are plotted in Fig. 8 (c). 


e 2 . 
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